First-principles calculation of magnetoelastic coefficients and magnetostriction 
in the spinel ferrites CoFe 2 4 and NiFe 2 4 
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We present calculations of magnetostriction constants for the spinel ferrites CoFe2C>4 and NiFe2 04 
using density functional theory within the GGA+(7 approach. Special emphasis is devoted to the 
influence of different possible cation distributions on the B site sublattice of the inverse spinel struc- 
ture on the calculated elastic and magnetoelastic constants. We show that the resulting symmetry- 
lowering has only a negligible effect on the elastic constants of both systems as well as on the mag- 
netoelastic response of NiFe2 04, whereas the magnetoelastic response of CoFe2 04 depends more 
strongly on the specific cation arrangement. In all cases our calculated magnetostriction constants 
are in good agreement with available experimental data. Our work thus paves the way for more 
detailed first-principles studies regarding the effect of stoichiometry and cation inversion on the 
magnetostrictive properties of spinel ferrites. 

PACS numbers: 75.80.+q, 71.15.Mb, 75.47.Lx 
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I. INTRODUCTION 

Magnetostriction describes the deformation of a 
ferro- or ferrimagnetic material during a magnetization 
process^— Thereby, one can distinguish between the 
spontaneous volume magnetostriction, which is indepen- 
dent of the magnetic field direction, and the so-called 
linear magnetostriction which characterizes the change 
of length along a certain direction that depends on the 
orientation of the applied magnetic field. The same mag- 
netoelastic interaction that causes magnetostriction also 
leads to changes in the magnetic anisotropy as function 
of an externally applied strain. 

Magnetostrictive materials are very important for 
applications as magnetic field sensors and magneto- 
mechanical actuators, where a large (and often also 
preferably linear) magnetic field response is essential. 8 
On the other hand magnetostriction also causes noise 
and frictional losses in magnetic transformer cores, so 
that in this context a minimization of magnetostriction 
is desirable. 

CoFe204 (CFO) is known to have one of the largest 
magnetostriction among magnetic materials that do not 
contain any resource-critical rare-earth elements^ It has 
thus recently come into focus for use in magnetostrictive- 
piezoelectric composites ) 10-12 where the goal is to achieve 
cross coupling between magnetic and dielectric degrees of 
freedom. Due to its insulating character and high mag- 
netic ordering temperature, CFO together with NiFe2 04 
(NFO) and other spinel ferrites is also a very attractive 
candidate for spintronics applications, in particular for 
spin-filtering tunnel barriersj 13 i 14 For many of these ap- 
plications, thin films of CFO and NFO are epitaxially 
grown on substrates with different lattice constants. The 
resulting substrate-induced strain can then lead to dis- 



tinctly different properties of the thin films compared to 
the corresponding bulk materials. 

In view of this, a good quantitative understanding of 
magnetoelastic properties of spinel ferrites, that provides 
a solid basis for the interpretation of experimental results 
and allows for further optimization of magnetostrictive 
properties, is highly desirable. In particular, the ability 
to accurately predict effects of cation off-stoichiometry 
or surface and interface effects can provide valuable in- 
sights into the fundamental mechanisms determining the 
observed properties. 

In previous work we have shown that first-principles 
calculations based on density-functional theory (DFT) 
provide a suitable description of the magnetoelastic prop- 
erties of spinel ferrites; 15 ' 16 thus demonstrating the feasi- 
bility of more detailed studies into strain-induced effects 
in thin film structures composed of CFO and NFO. Here 
we extend our previous study, in order to provide a more 
comprehensive picture of the magnetoelastic response of 
CFO and NFO, in particular including first-principles 
calculations of the complete set of cubic magnetoelastic 
and magnetostrictive coefficients. Most importantly, we 
investigate the influence of different possible cation dis- 
tributions on the spinel B site sublattice on the magne- 
toelastic response of these materials. The purpose of the 
present work is to provide a first-principles based descrip- 
tion of magnetoelastic coupling in spinel ferrites that can 
be used as basis for further studies of the effect of cation 
substitution or off-stoichiometry on the magnetostrictive 
properties of this important class of materials. 

This paper is organized as follows. In Sec. Ill Al the 
spinel crystal structure is discussed, with special empha- 
sis on cation inversion and different possible cation ar- 
rangements on the B site sublattice. A general overview 
of magnetoelastic theory in cubic and tetragonal crystals 
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FIG. 1: (Color online) The spinel structure consists of an fee 
network of oxygen anions (red) with cations occupying differ- 
ent interstitial sites of the fee lattice, resulting in tetrahedrally 
coordinated A sites (purple) and octahedrally coordinated B 
sites (brown). Picture has been generated using VESTA. 17 



is given in Sec. Ill Bl Sec. Ill Cl describes how we determine 
all elastic and magneto-elastic coefficients from total en- 
ergy electronic structure calculations, while Sec. lIIDl pro- 
vides some more technical details of our calculations. Our 
results for CFO and NFO are presented in Sec. IIIII and 
our main conclusions are summarized in Sec. II VI 



II. 



THEORETICAL BACKGROUND AND 
COMPUTATIONAL METHOD 

Inverse spinel structure and different cation 
distributions 



for short-range cation order on the B sites have been re- 
ported recently for the case of NFO, both in bulk single 
crystals as well as in thin filmsi 19 ' 20 

In the present work we represent the inverse spinel 
structure within a tetragonal unit cell containing 4 for- 
mula units (see also Ref. |21| ) using lattice vectors a*i = 
(a/2, -a/2,0), a 2 = (a/2, a/2,0), and a 3 = (0,0, c), so 
that c/a = 1 corresponds to the unstrained, nominally 
cubic case. By distributing equal amounts of Co (respec- 
tively Ni) and Fe on the 8 B sites within this unit cell, 
70 cation arrangements belonging to 8 different space- 
groups can be generated. In the following we consider 
only the three high-symmetry arrangements shown in 
Fig. [2] (a)-(c), plus one additional low-energy configu- 
ration for CFO, corresponding to 75% inversion, shown 
in Fig. [5] (d). The specific cation arrangements shown in 
Fig. [5] in combination with the periodic boundary con- 
ditions corresponding to the tetragonal lattice vectors 
reduce the space group symmetries to P4i22 (No. 91), 
Imma (No. 74), and P4m2 (No. 115) for the fully inverse 
configurations, and to PI (No. 1) for the case with 75 % 
inversion. As we have previously shown^ both P4i22 
and Imma correspond to low energy configurations for 
the fully inverse case, with P4i22 slightly lower in energy 
than Imma for both CFO and NFO, whereas the P4m2 
configuration is energetically much less favorable. The 
PI structure represents a low energy configuration for the 
case A = 0.75i2i We also note that the P4i22 configura- 
tion corresponds to the local structure suggested for the 
experimentally observed short-range order in NFO ) 19 i 20 
whereas the Imma configuration is equivalent to the one 
used in our previous study of magneto-elastic effects in 
CFO and NFO^i£ 



Both CFO and NFO crystallize in the cubic spinel 
structure (see Fig. [T]), which belongs to space group 
Fd3m (No. 227). The spinel structure contains two in- 
equivalent cation sites, a tetrahedrally coordinated A site 
and an octahedrally coordinated B site. In the normal 
spinel structure each of these sites is occupied by a par- 
ticular cation species (e.g. divalent Mn 2+ on the A site 
and trivalent Fe 3+ on the B site in the case of MnFe204). 
However, in the inverse spinel structure, the more abun- 
dant cation species (here: Fe 3+ ) occupies all A sites and 
50 % of the B sites, with the remaining 50 % of B sites 
occupied by the less abundant cation species (here: Co 2+ 
or Ni 2+ ). In practice, intermediate cases can also occur, 
characterized by an inversion parameter A, ranging from 
A = for the normal spinel structure to A = 1 for com- 
plete inversion. 

Both CFO and NFO are experimentally found to be 
inverse spinels, with Awl for NFO but only incom- 
plete inversion for CFO (with A between 0.76 — 0.93, de- 
pending strongly on sample preparation conditions) 
Both materials are generally found to be perfectly cu- 
bic, with a random distribution of divalent and trivalent 
cations over the B site sublattice. However, indications 



B. Magnetoelastic theory 

Within the phenomenological theory of magnetoelas- 
ticity, the magnetoelastic energy density / = E/V is 
expressed in terms of the direction cosines of the magne- 
tization vector, cti (i = x,y,z), and the components of 
the strain tensor ey, relative to a suitably chosen (non- 
magnetic) reference stated— This energy density can be 
divided into a purely elastic term, f e \, and a magnetoe- 
lastic coupling term, / mc , which is usually taken as linear 
in the strain components. For a cubic crystal these terms 
have the following form^ 

ycubic =-Cn(e xx _|_ £ 2^ _|_ £ _.j _|_ 2Ca{e 2 xy + e'y Z + e 2 zx ) 

~r"Cl2 {^-yy^-zz "t" &xx^zz "T" ^xx^-yy) j 

(1) 

and 

+B 1 (a 2 x £xx H~ &y£yy ~l~ ^2^22) (2) 

-\- i 2B2^Ct X (XyE X y -\~ CXyOi Z Sy Z ~\~ OL Z Oi X £ Z X ) 7 



FIG. 2: (Color online) Cation distribution of Fe (brown) and Co (Ni) (blue) on the B sites of the spinel structure for the 
different configurations used in our calculations. Note that only the B sublattice is shown. From left to right the depicted 
structures correspond to spacegroups (a) P4i22 (No. 91), (b) Imma (No. 74), and (c) P4m2 (No. 115). Figure (d) on the 
right displays the CFO low-energy solution with incomplete degree of inversion, A = 0.75 corresponding to spacegroup PI (No. 
1).— Pictures have been generated using VESTA.— 



where Cn, C\%, and C44 are elastic and Bq, B\, and B2 
are magnetoelastic coupling constants. 

The relative length change along an arbitrary (measur- 
ing) direction with direction cosines /3j is given by: 



(3) 



where the strain components depend on the magnetiza- 
tion directions. These equilibrium strains as function of 
the magnetization direction can be found by minimizing 
the sum of the two energy expressions (Q]) and © with 
respect to all strain components. This results in: 



Al 

T 



cubii 



=A" + -A 100 (a 2 J 2 x + al(3 2 y + a 2 J 2 z - y 3 ) 



+3Ain (a x a y p x fi y + a y a z /3 y fi z + a x a z /3 x {3 z ) ■ 

(4) 

Here, A" = -(B + B 1 /3)/{C 11 + 2C 12 ) describes a pure 
volume magnetostriction that is independent of the mag- 
netization direction (this term is sometimes omitted from 
the above formula and is of no concern in the present 
work). The widely used magnetostriction constants of a 
cubic crystal are given by: 
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Am — — 
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(5) 



(6) 



These two coefficients measure the fractional length 
change along the [100] [fi x = l,(3 y = (3 Z = 0) and [111] 
(/3j = l/v3) directions, when the sample is magnetized 
to saturation along the [100] (a x = 1, a y = a z = 0) and 
[111] {on — l/\/3) directions, relative to an ideal demag- 
netized reference state which is defined by (oif) = V3 
and (pnot.j) — 0. In a polycrystalline sample one can only 



measure a direction average over both A100 and Am given 
by* 



Ac 



2^ 3^ 
g A100 + g Am • 



(7) 



As noticed in Sec. Ill Al the cation arrangements used 
to describe the inverse spinel structure within our cal- 
culations lower the cubic symmetry of the ideal spinel 
structure to tetragonal (P4i22 and P4m2), orthorhom- 
bic (Imma), or even triclinic (PI). A full first-principles 
description of magnetoelastic effects within these lower 
symmetries would require the calculation of 6 (9, 21) 
different elastic and 7 (12, 36) magnetoelastic coupling 
constants for the mentioned tetragonal (orthorhombic, 
triclinic) spacegroups, respectively^ Due to the resulting 
large computational effort, and considering the fact that 
experimentally both CFO and NFO are found to be cu- 
bic, we do not attempt such a full determination of all 
elastic and magnetoelastic coefficients within the lower 
symmetries, and instead evaluate our results using the 
relations for the cubic case described above (i.e., similar 
to our previous work in Refs. [TJ| and ITfjl). To estimate 
the degree to which the lower symmetry affects our cal- 
culated coefficients, we also compare some of our data 
to the correct formulas corresponding to the lower sym- 
metry. For simplicity we hereby restrict ourselves to the 
tetragonal case. The required equations are presented in 
the following. 

Within the lower tetragonal symmetry there are six 
independent elastic and seven different magnetoelastic 
coupling constants, in contrast to the three elastic and 
three magnetoelastic coefficients in the cubic case^ The 
resulting expressions for f c \ and / me then read£ 



ftct _ 1 ( 2 
/el - 2 C H \ £ xx + e 



2 ) 

mi) 



1 



; c 33£ 2 



Cl2£xx£yy ^13 (_€ X x H~ &yy) &zz 
„2 



(8) 



2c 44 (e 2 yz 



2C66£ 



xy ' 



with aj denoting the six different tetragonal elastic con- 
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stants, and 

/me =bll (£xx + £yy) + b\ 2 £ zz 

+b 21 (al - V 3 ) (Exx + Eyy) + b 22 (a 2 z - V 3 ) 
+\ b z (a 2 x - al) (e xx - £ yy ) + b 3 a x a y e xy 
+& 4 (a x a z £ xz + a y a z £ yz ) , 



(9) 



with the various 6's denoting the seven different tetrago- 
nal magnetoelastic coupling constants. The correspond- 
ing cubic expressions (JTJ) and ^ can then be obtained 
from ([8]) and ([9]) with the additional symmetry con- 
straints: en = c 33 = C11, C12 = ci 3 = C12, c 44 = c 66 = 
C44, foil = b 12 = B + V3-B1, b 22 = -2b 21 = 63 = Si, 
and 63 = 64 = B 2 . 



C. Determination of elastic and magnetoelastic 
constants 

In order to determine the (cubic) elastic constants 
for CFO and NFO, we first perform a full structural 
relaxation of both systems. Similar to our previous 
investigations) 15 ^ 6 ' 21 we thereby constrain the lattice 
vectors to "cubic" symmetry (c/a = 1) and fix the in- 
ternal coordinates of the A and B cations to ideal values 
corresponding to the cubic spinel structure, i.e., we only 
allow for an optimization of the total volume and the 
oxygen positions. We then determine the three indepen- 
dent cubic elastic constants Cu, C12, and C44, and the 
two cubic magnetoelastic coupling constants B\ and B 2 
by distorting the equilibrium crystal structure in three 
different ways: i) isotropic volume expansion, ii) con- 
straining two of the three lattice dimensions and relax- 
ing the third ("epitaxial strain"), and iii) by applying a 
volume-conserving shear strain. 

i) Isotropic volume expansion. The dependence of the 
total energy E to t on the unit cell volume V provides the 
bulk modulus B, which is defined as 



B = V 



d 2 E toi 
dV 2 



(10) 



(V=V ) 



with Vq being the equilibrium volume. According to 
Eq. ((T|) the bulk modulus B of a cubic crystal can be 
expressed in terms of the elastic moduli Cu and C12 as 
follows: 



B= X - (C u + 2C 12 ) 



(11) 



ii) Epitaxial strain. We follow the approach of Ref. |15| to 
obtain a second independent elastic constant by applying 
epitaxial strain, i.e., we constrain the "in-plane" lattice 
constant to values ranging from —4% to +4% relative to 
the theoretical equilibrium lattice constant ao, and we 
relax the "out-of-plane" lattice constant and all internal 



structural parameters of the oxygen anions. The rela- 
tion between the relaxed out-of-plane strain e±_ and the 
fixed in-plane strain en then defines the so-called two- 
dimensional Poisson ratio v 2 r>. It follows from Eq. ([1]) 
that for a cubic system v 2 £> is given as: 



V 2 D 



£-L _ 2 ^i2 
£|| Cu 



(12) 



The elastic moduli Cu and C12 can then be obtained 
from Eqs. (ITT|) and (fT2"|) using the bulk modulus and two- 
dimensional Poisson ratio calculated from DFT. 

For the cation arrangements with tetragonal, or- 
thorhombic, or triclinic symmetry depicted in Fig. [5] the 
ratio £_i_/ e || can be different for different orientations of 
"out-of-plane" and "in-plane" directions relative to the 
crystal axes. To quantify the resulting difference we per- 
form calculations for two symmetry-inequivalent orien- 
tations of the applied strain em. In particular we ap- 
ply the epitaxial constraint first within the xy plane 
{£\\ = £n = £yy and e± = £ zz ) and then also within 
the yz plane (e\\ = £ yy = £ zz and e± = e xx ). Using 
the tetragonal energy expressions of Eqs. ([5]) and ^ to- 
gether with the definition of v 2 d in Eq. (fT2|) one obtains 

Vtfp = 2ci 3 /c 3 3 and i/^n = ( c i2 + ci 3 )/cn for these two 
cases. The difference between these two values for v 2 u 
thus gives a measure for the difference between cu and 
C33 as well as between C12 and C13. 

To obtain the magnetoelastic coupling coefficient B\ 
we monitor the total energy differences for different orien- 
tations of the magnetization as a function of the applied 
in-plane constraint en and relaxed out-of-plane strain 
£± = —v 2 r,£\\. Using the cubic expression ([2]) for / mG one 
can see that the strain dependence of the energy density 
for all in-plane orientations of the magnetization is given 
by B\ ■ e\\ , whereas the strain dependence for out-of-plane 
orientation is given by —B± ■ v 2 r> ■ em. The strain depen- 
dence of the total energy difference between out-of-plane 
versus in-plane orientation of the magnetization is thus 
given byj^ 3 . 



AE/V=-(p2D + l)B 1 e ]{ . 



(13) 



The coefficient B\ can therefore be obtained from the 
calculated strain-dependent magnetic anisotropy energies 
(MAEs) and the previously determined two-dimensional 
Poisson ratio v 2 r> ■ While B\ is not directly accessible by 
experimental investigations, it is related to the magne- 
tostriction constant A100 via Eq. ([5]). 

In the tetragonal case the monitored strain dependence 
of the total energy difference between out-of-plane ver- 
sus in-plane directions of the magnetization will depend 
on the orientation of "out-of-plane" and "in-plane" direc- 
tions with respect to the tetragonal crystal axes. For the 
epitaxial constraint applied within the xy plane (i.e. £|| = 

£xx — ^yyi 

leading to a Poisson ratio v^f — 2ci3/c 33 ) 
and using the tetragonal energy density (Eqs. ([5]) and 
(|9"|)), the following expression for the strain dependence 
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of the total energy difference between in-plane and out- 
of-plane magnetization can be obtained: 

(AE)^/V = (26a! - ^ ) 6 22 )£|| , (14) 
which is valid for all in-plane orientations of the mag- 



netization. In contrast, for the epitaxial constraint ap- 
plied within the yz plane (i.e. en = e yy — e zz , leading 

to a Poisson ratio v^j^ = {c.\% + ci^)/cn) the resulting 
(AE)^^ /V depends on the specific in-plane direction 
and is given by: 



u (,,(v z ) 



(AE) (yz) /V = 



■f 1) - (6ai + &22 



,(v z 



V 2D b 21 



-<>3{y 2D 



..(yz 

y 2D 



hi. 



for (A£)^« 


— £100 


— Eqqi 


for (AE) ( y^ 


= £a 00 


— Eqiq 


for (A£)^ 2 ' 


= £100 


~ £011/011 



(15) 



iii) Volume- conserving shear strain. The third cubic 
elastic modulus C44 is calculated according to Mchl, 24 
by applying a volume-conserving monoclinic shear strain 
in the xy plane (e\\ = e xy , e± = e zz = e|/(l - e|), 

£xx = £yy = £yz — e zx = 0). The resulting change in 
total energy can then be written as: 

E(±e\\) = WC u e\ + 0[e{] , (16) 

which allows for a straight-forward determination of C44. 

For the cation arrangements with tetragonal, or- 
thorhombic, or triclinic symmetry depicted in Fig. [2] 
different shear planes (s xy , e yz , e zx ) are connected to 
different elastic moduli cu. Using the tetragonal en- 
ergy expressions of Eqs. © and ([5} together with the 
volume-conserving monoclinic strain in the xy plane de- 
scribed above, one notices the connection of e xy and cqq. 
However, choosing a volume-conserving monoclinic strain 
in the yz plane (£|| = e yz , e± = e xx = e|/(l - ejj), 

£yy = Ezz = £xy = £zx = 0) yields directly C44, allowing 
for a comparison with cqq. 

Similar to the first magnetoelastic coupling constant 
B%, the second coefficient B% is determined by monitor- 
ing the total energy differences between different orien- 
tations of the magnetization as a function of the applied 
strain en. Depending on whether the shear strain e\\ is 
applied within the xy or yz plane, we consider the fol- 
lowing energy differences: 

(AE)^ = E no - £100/010 = £100/010 _ £110 (17) 

(AEp yz ' = E ii — £010/001 = £010/001 — £011 • (18) 

In all cases, the strain-dependence of these total energy 
differences can be written as: 

AE/V = B 2 e l{ . (19) 

Thus, the strain dependence of these energy differences 
is governed by the magnetoelastic coupling constant £2 > 
which can be determined from the calculated A£/V(e||). 

Similar to £1 , the magnetoelastic coupling constant £2 
is also not directly accessible by experiment, but it is re- 
lated to the magnetostriction constant Am via Eq. ([5]). 



Once the magnetostriction constants A100 and Am are 
obtained, the average magnetostriction constant Ag, suit- 
able for polycrystalline samples, can be calculated from 
Eq. 0. 



D. Other computational details 

All calculations presented in this work are performed 
using the projector-augmented wave (PAW) method^ 
implemented in the Vienna ah initio simulation package 
(VASP 4.6)^^ Standard PAW potentials supplied with 
VASP were used in the calculations, contributing nine va- 
lence electrons per Co (4s 2 3d 7 ), 16 valence electrons per 
Ni (3p 6 4s 2 3d 8 ), 14 valence electrons per Fe (3p 6 4s 2 3d 6 ), 
and 6 valence electrons per O (2s 2 2p 4 ). 

The generalized gradient approximation according to 
Perdew, Burke, and Ernzerhof (PBE)2£i is used in com- 
bination with the Hubbard "+J7" correction, 31 where 
U=3 eV and J—0 eV is applied to the d states on all 
transition metal cations. We have shown in Refs. [HI, 
[HI, and [21] that this gives a realistic description of the 
electronic structure of CFO and NFO and leads to re- 
sults which are in good overall agreement with available 
experimental data. 

All structural relaxations are performed within a 
scalar-relativistic approximation, whereas spin-orbit cou- 
pling is included for the calculation of the MAEs. A plane 
wave energy cutoff of 500 eV is used, and the Brillouin 
zone is sampled using a T-centered 5x5x3 fc-point grid 
both for the structural optimization and for all total en- 
ergy calculations. We have verified that all quantities of 
interest, in particular the magnetic anisotropy energies, 
are well converged for this /c-point grid and planewave 
energy cutoff. 
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III. RESULTS AND DISCUSSION 



A. Structural properties 



The equilibrium lattice constants do, bulk moduli B, 
two-dimensional Poisson ratios V2D, the resulting elastic 
constants C\\ and C12, as well as C44 obtained for the dif- 
ferent cation arrangements for both CFO and NFO are 
given in Tab. HI One notices that the calculated lattice 
constants for the two low-energy configurations Imma 
and P4i22 are very similar to each other, and that the 
ones for the higher energy P4m2 configuration and for 
the case with 75 % inversion for CFO are slightly larger 
than that (by less than 0.2%). This increase in lattice 
constant is mirrored by a corresponding decrease in the 
bulk modulus (by about 3%). Overall, the variation 
of both bulk modulus and equilibrium lattice constant 
between different cation distributions is much smaller 
than the slight under- and overestimation of these quan- 
tities with respect to the experimental value, which is 
within the usual limits of the PBE+C/ approach (see also 
Ref.QJ). 

It can also be seen that the difference in the two- 
dimensional Poisson ratios obtained for two different ori- 
entations of e± is rather small and of similar magnitude 
as the differences between the various cation arrange- 
ments. This indicates that the symmetry-lowering due to 
the different cation arrangements has only a small effect 
on the elastic properties, which can still to a good ap- 
proximation be described by cubic elastic constants C\\ 
and Ci2- 

Applying the volume-conserving monoclinic strain as 
described in Sec. IIIBI yields the remaining elastic mod- 
ulus C44 which is in very good agreement with the ex- 
perimental values for both CFO and NFO. To evaluate 
the influence of different orientations of e±_ on C44 we 
applied e± ~ e xy and e±_ — e yz with the respective £11 
to the low energy orthorhombic Imma symmetry. The 
difference in the obtained C44 is slightly larger compared 
to the difference in the C\\ and C12, but still within the 
typical uncertainties of first-principles methods. 

Overall it appears that while the agreement between 
the calculated and experimental lattice constants and 
clastic moduli is quite good and within the typical un- 
certainties of state-of-the-art first-principles methods, 
the uncertainties resulting from the symmetry-lowering 
cation arrangements are significantly smaller than that. 
Therefore, the elastic properties of the various cation ar- 
rangements of lower symmetry can be well described by 
cubic elastic constants. 
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FIG. 3: Total energy difference AE per two formula units 
(f.u.) of NFO as function of the epitaxial constraint eii for 
different cation arrangements. The left (right) panels corre- 
spond to the case with e± — e zz (e± = e xx )- The panels from 
top to bottom refer to symmetries P4i22, Imma, and P4m2, 
respectively. In case of e±_ = e Z z [s± = e xx ) the depicted en- 
ergy difference AE is taken with respect to the [001] ([100]) 
direction, with the symbols denoting A [100] ([010]), T [010] 
([001]), < [110] ([Oil]), and ► [110] ([011]), respectively. 



B. Magnetoelastic properties 

1. NFO 

Next we focus on the magnetoelastic coupling in NFO. 
The calculated MAEs necessary to determine the mag- 
netoelastic coupling constant B\ are depicted in Fig. [3] 
As described in Sec. Ill CI these MAEs are defined here 
as the energy differences for various orientations of the 
magnetization with respect to the magnetization direc- 
tion perpendicular to the applied strain plane, i.e., [001] 
for £11 = e xx = £ yy and [100] for £11 = e yy = e zz . Accord- 
ing to Eq. ([T3")l the slope of the curves given in Fig. [3J is 
directly related to the magnetoelastic coupling constants 
B\. At first sight, the slopes of all curves in all panels 
are very similar and negative, thus leading to a positive 
B\ (the range of the y axes is the same in all panels to 
allow for a direct inspection of slope differences) . 

In the tetragonal symmetries (P4i22 and P4m2) all 
curves fall on top of each other for e± = e zz (left pan- 
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TABLE I: Optimized equilibrium lattice constant ao, bulk modulus B, two-dimensional Poisson ratio U2D, and elastic moduli 
Cn, Ci2, and C44 for CFO and NFO, obtained for different cation arrangements and strain orientations (ex = £ zz = z 
and ex = e xx = a;) in comparison to experimental data. The experimental V2D has been evaluated from Eq. (|12[1 using 
the experimental elastic constants. PI in case of CFO refers to the low-energy solution with incomplete degree of inversion, 
A = 0.75.— 



CFO 


ao 


B 




V2D 


On 


O12 


O44 


(A) 


(GPa) 




(GPa) 


(GPa) 


(GPa) 


Imma 


8.463 


172.3 


z 

X 


1.132 
1.147 


242.5 
240.8 


137.3 
138.1 


94.9 
83.2 


P4i22 


8.464 


170.8 


X 


1.129 
1.147 


240.7 
238.7 


135.9 
136.9 


84.7 


P4m2 


8.473 


168.0 


X 


1.132 
1.128 


236.4 
236.8 


133.8 
133.6 


92.3 


PI 


8.477 


167.8 


Z 


1.155 


233.6 


134.9 


87.7 


X 


1.146 


234.6 


134.4 




Exp. (Ref. 32) 


8.392 


185.7 




1.167 


257.1 


150.0 


85.3 


NFO 


ao 

(A) 


B 
(GPa) 


e± 


V2D 


Cn 
(GPa) 


C12 
(GPa) 


C44 
(GPa) 


Imma 


8.426 


177.1 


z 

X 


1.106 
1.115 


252.2 
251.2 


139.5 
140.0 


93.2 
87.6 


P4i22 


8.428 


175.4 


z 


1.116 


248.7 


138.8 


87.4 


X 


1.116 


248.7 


138.8 




P4m2 


8.435 


173.3 


z 

X 


1.116 
1.102 


245.7 
247.3 


137.1 
136.3 


91.0 


Exp. (Ref. 32) 


8.339 


198.2 




1.177 


273.1 


160.7 


82.3 



TABLE II: Magnetoelastic coupling constants (Pi, B2) and 
magnetostriction constants (A100, Am, As) for NFO using 
different cation arrangements and strain planes according to 
e± — Ezz — z (ex = £xx = x) in comparison with available 
experimental data. The average magnetostriction constant 
As has been obtained using Eq. ©. 







Bi 


A100 


B 2 


Am 


As 






(MPa) 


(xlO" 6 ) 


(MPa) 


(xlO 


- 6) 


P4i22 


z 

X 


6.6 
6.4 


-40.1 
-38.6 


2.5 


-9.7 


-21.9 


Imma 


z 


6.1 


-35.9 


0.9 


-3.4 


-16.4 


X 


6.7 


-40.3 


1.9 


-7.3 


-20.5 


PAm2 


z 

X 


6.5 
6.9 


-40.0 
-41.3 


1.4 


-5.3 


-19.2 




Ref. 33° 




-36.0 




-4.0 


-16.8 


Exp. 


Ref. U b 




-50.9 




-23.8 


-34.6 




Ref. 35 c 




-43.0 




-20.1 


-29.3 



"Single crystals with Nio.sFe2.2O4 composition. 
''Single crystals of NiFe2C>4. 
c Single crystals of NiFe204. 



els), whereas there is a small offset between the curves 
in all other cases, due to the lower symmetry. In the 
even lower Imma symmetry this offset is also present for 
£1 = s Z z- Nevertheless, the variation with strain is very 
similar in all cases, and the values for B%, obtained by 
averaging over all curves corresponding to the same sym- 
metry and strain orientation, are given in Tab. [Til These 
values range from 6.1 MPa to 6.9 MPa, depending on the 
specific cation arrangement and strain orientation. Due 



to these rather small variations, we can conclude that the 
magnetostrictive response in NFO can to a good approx- 
imation be described as cubic. 

Together with the respective elastic constants from 
Tab. U the magnetostriction constants A100 can be ob- 
tained via Eq. <j5j> , and are also listed in Table HU It 
can be seen that there is only a weak influence of ei- 
ther cation arrangement or different strain planes on the 
NFO magnetostriction constant A100, which ranges from 
-35.9 x KT 6 to -41.3 x 1CT 6 . This agrees perfectly 
with experimental data ranging from —36.0 x 10~ 6 to 
-50.9 x 10~ 6 . 

The calculated strain-dependent MAEs necessary for 
the determination of B2 are shown in Fig. |U The dif- 
ferent curves are adjusted to match at e\\ = in order 
to remove the corresponding offset which is irrelevant for 
the present work. The MAEs are chosen according to 
Eqs. (jTTJ) and (fT5]) as energy differences between different 
in-plane orientations of the magnetization with respect 
to the applied shear strain eii . According to Eq. (IT9l) the 
slope of the curves given in Fig. @] is directly related to 
the magnetoelastic coupling constant B 2 - From Fig. 0] 
it can be seen that the slopes of these curves are pos- 
itive, corresponding to positive B 2 . There are slightly 
stronger nonlinearities in the curves in each of the panels 
compared to Fig. [31 as well as a stronger influence of the 
explicit cation arrangements. The resulting magnetoelas- 
tic coupling constants B 2 are listed in Tab. Inland range 
from 0.9 MPa to 2.5 MPa, leading to magnetostriction 
constants Am ranging from —3.4 x 10~ 6 to —9.7 x 10~ 6 . 
These values are compatible with the lower experimen- 
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[%] 




0.05 ^ 



FIG. 4: Total energy difference AE per two formula units (f.u.) of NFO as function of shear strain for different cation 
arrangements. The panels from left to right refer to symmetries P4i22 (en = Exy), Imma (eii = Exy), Imma (en = E yz ), and 
P4m2 (ey = e xy )_, respectively. The depicted energy differences AE correspond to [110]-[100] (A), [110]-[010] (T), [100]-[ll0] 
(►), and [010]-[110] (■<), respectively, for en = e xy and equivalent directions for en = e yz . 



tal values, which themselves range from —4.0 x 10 6 to 
-23.8 x 10~ 6 . The last column in Tab. HQ also lists the 
averaged As suitable for polycrystalline materials using 
Eq. 0. 

Overall, the different cation arrangements and strain 
planes have only a rather weak influence on the calculated 
magnetostriction constants of NFO, which agree very 
well with the range of reported experimental data. We 
can therefore confirm our earlier finding^ that BFT+U 
methods are suitable for a quantitative description of 
magnetoelastic properties in this material. Moreover, 
although the symmetries of the investigated cation ar- 
rangements are not cubic, the magnetostrictive proper- 
ties of NFO are very well described within the cubic the- 
ory. 



2. CFO 

Now we turn to our results for CFO. The calculated 
MAEs for the determination of the magnetoelastic cou- 
pling constant B\ are depicted in Fig. [SJ At first sight, 
one notices again that all slopes are negative, leading to 
a positive magnetoelastic coupling constant B\. How- 
ever, in contrast to NFO, the values are now much larger 
and also depend more strongly on the specific cation ar- 
rangement and orientation of the strain plane. In all 
cases except for the case of P4m2 with eii ~ e xx , we 
again obtain an offset between the different curves, which 
is due to the lower symmetry of the specific cation dis- 
tribution. The differences in slopes observable between 
the various curves in the left panel of P4i22 symmetry 
are due to the fact that in this case the system adopts 
an orbitally-ordered ground state with symmetry lower 
than that of the underlying crystal structure. Strongest 
deviations from linearity are observed in the low-energy 
solution with symmetry PI belonging to incomplete in- 
version A = 0.75. 

The determined magnetoelastic coupling constants B\ 
are given in Tab. QUI ranging from 18.9 MPa to 42.0 MPa. 
The largest influence of the strain plane orientation is ob- 
served for P4i22 symmetry. Overall, the specific cation 



TABLE III: Magnetoelastic coupling constants (Pi, P2) and 
magnetostriction constants (A100, Am, As) for CFO using 
different cation arrangements and strain planes according to 
e± = E zz = z (e± = e X x = x) in comparison with available 
experimental data. The average magnetostriction constant 
As has been obtained using Eq. (J7J). 







Bi 


A100 


B 2 


Am 


As 




£± 


(MPa) 


(xlO" 6 ) 


(MPa) 


(xl0~ 


- B) 


P4i22 


z 


18.9 


-120.1 


-8.4 


32.9 - 


-28.3 


X 


32.8 


-215.0 








Imma 


z 


39.7 


-251.7 


-11.6 


40.9 - 


-76.1 


X 


29.2 


-189.7 


-12.2 


48.8 




Pim2 


z 


42.0 


-272.7 


-14.6 


52.7 - 


-77.5 


X 


30.9 


-199.4 








PI 


z 


29.1 


-196.3 


-7.6 


28.8 - 


-61.2 


X 


24.2 


-160.9 










Ref . 36 a 




-225.0 








Exp. 


Ref. 33 b 




-250.0 










Ref. 33 c 




-590.0 




120.0 - 


164.0 



"Polycrystalline CoFe2 04. 

''Single crystals with C01.1Fe1.9O4 composition. 
c Single crystals with Coo.8Fe2.2O4 composition. 



arrangement has a much larger influence on the obtained 
magnetoelastic coupling constants in CFO compared to 
NFO. However, we note that even though there are pro- 
nounced differences between the two different strain ori- 
entations (left and right panels in Fig. [5]) for the same 
cation arrangements, the strain dependence of the vari- 
ous calculated energy differences for the same strain ori- 
entation (different curves within each panel) are very 
similar in each case. From expression (|15l) for tetrago- 
nal symmetry, we can therefore empirically observe that 
the following approximate relationship holds between the 
various magnetoelastic coefficients: 



1 



b-i(v 2D + 1) ~ b 2 i + b 22 - v w bix 



(20) 



However, since the slopes in the left and right panels of 
Fig. differ, the stronger condition 6 3 = b 22 = — 26 2 i: 
which would be valid within cubic symmetry, is not ful- 
filled in CFO. The deviation from cubic symmetry caused 
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FIG. 5: Total energy difference AE per two formula units 
(f.u.) of CFO as function of the epitaxial constraint sii for 
different cation arrangements. The left (right) panels corre- 
spond to the case with e± = e zz (e± = e xx ). The panels from 
top to bottom refer to symmetries P4i22, Imma, and P4m2, 
and PI (low-energy solution for cation inversion A = 0.752iY 
In case of e± = e zz (s± = e xx ) the depicted energy difference 
AE is taken with respect to the [001] ([100]) direction, with 
the symbols denoting A [100] ([010]), T [010] ([001]), <4 [110] 
([011]), and ► [110] ([011]), respectively. 



The strain-dependent MAEs necessary for the determi- 
nation of B2 are shown in Fig. [5J analogous to the NFO 
case. Most strikingly, and in contrast to NFO, the corre- 
sponding slope is negative, thus leading to a negative Bi 
in CFO. The spread in slopes in each of the panels is com- 
parable to NFO. While we obtain quite similar values for 
Imma and P4m2 symmetry (middle three panels), and 
also for P4i22 and PI, the latter two symmetries lead to 
somewhat smaller values for Bi than the former. 

Overall, B 2 ranges from -8.4 MPa to -14.6 MPa for 
the symmetries corresponding to complete cation inver- 
sion, and -7.6 MPa for the case with A = 0.75 (PI). The 
resulting magnetostriction constants Am of CFO range 
from 28.8 x 10 -6 to 52.7 x 1CT 6 , respectively. These val- 
ues are lower than the (to the best of our knowledge only 
available) value of 120 x 10~ 6 reported experimentally. 

In view of the relatively strong dependence on the spe- 
cific cation arrangement, no particular trend is apparent 
on how the magnetostriction constants change with re- 
duced cation inversion (PI structure compared to the 
other cases with full inversion). Taking a closer look at 
the individual magnetostriction constants Aioo for all in- 
vestigated cation arrangements and strain planes, one 
can notice that the largest magnetostriction occurs for 
cases where the cation species are arranged in alter- 
nating planes parallel to the applied strain plane, e.g., 
£± = £zz — z for Imma and P4m2 symmetry (see Fig. [2]). 
Furthermore, if one compares the two different strain ori- 
entations for P4i22 symmetry, the magnetostriction is 
larger for e± = e xx = x, where the strain plane contains 
chains of B site cations with two equal cations next to 
each other in each chain. The magnetostriction value 
for the strain plane containing alternating cation chains 
within P4i22 symmetry is the smallest observed here. 
However, at present it is unclear whether these correla- 
tions between cation arrangement and Aioo are mostly 
coincidental, or whether they indeed indicate a deeper 
relationship between these two properties. In any case 
our results give clear evidence that a fully quantitative 
model of anisotropy and magnetostriction in CFO needs 
to include crystal- or ligand-field effects that go beyond 
the immediate nearest neighbor shell of the Co 2+ cation. 



by the specific cation arrangements, is therefore more 
strongly manifested in the magnetoelastic response of 
CFO compared to NFO. Nevertheless, the approximate 
relation Eq. (|20|) indicates that some residue of the ap- 
proximate structural cubic symmetry is still present also 
in the case of CFO. 

The magnetostriction constants of CFO can now be 
obtained via Eq. ([5]) and using the elastic constants in 
Table Q] The resulting values are listed in Table IIIII and 
range from -120.1 x 10~ 6 to -272.7 x 10 -6 . This agrees 
well with the lower range of available experimental data, 
which itself varies between —225 x 10 -6 and —590 x 10~ 6 . 



The effect of different distributions of Co 2+ and Fe 3+ 
cations on the B sites surrounding a specific Co B site 
has been taken into account in the theory of magnetic 
anisotropy for CFO by Tachikij^I and is also discussed 
by Slonczewskij^ It was shown that the correspond- 
ing crystal-field component can have a strong effect on 
the resulting cubic magnetic anisotropy constants. The 
noticeable dependence of our calculated magnetoelastic 
coupling constants on the specific cation arrangement in 
CFO indicates that this crystal-field component is indeed 
quite strong and needs to be taken into account within 
a quantitative theory of anisotropy and magnetostriction 
in spinel ferrites. 
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FIG. 6: Total energy difference AE per two formula units (f.u.) of CFO as function of shear strain for different cation 
arrangements. The panels from left to right refer to symmetries P4i22 (su = e xy ), Imma (em = e xy ), Imma (e« = e yz ), P4m2 
(e|| = £xy), and PI (e|| = e xy , low-energy solution for cation inversion A = 0.75^), respectively. The depicted energy differences 
AE correspond to [110]-[100] (▲), [110]-[010] (T), [100]-[110] (►), and [010]-[110] (<), respectively, for ey = e xy and equivalent 
directions for ey = e yz . 



IV. SUMMARY AND CONCLUSIONS 

In summary, we have presented a detailed first- 
principles study of elastic and magnetoelastic properties 
of the inverse spinel ferrites NFO and CFO. We have 
calculated all cubic elastic and magnetoelastic constants 
from a variety of distorted crystal structures. Thereby, 
we have considered different possible cation arrangements 
to represent the inverse spinel structure, and in the case 
of CFO we also considered a cation distribution corre- 
sponding to incomplete inversion with A = 0.75. The 
magnetoelastic coefficients are obtained from the strain 
dependence of the MAEs for two different deformations 
of the crystal structure. 

Even though the symmetry of the considered cation 
arrangements is lower than cubic, our results show that 
the elastic response of both NFO and CFO can to a good 
approximation be described using cubic elastic constants. 
Since the elastic constants are mainly determined by the 
strength of the chemical bonding, this indicates that Co, 
Ni, and Fe all form bonds of similar strength with the 
surrounding atoms. 

Similarly, the magnetoelastic response of NFO can also 
to a good approximation be described using the cubic ex- 
pression for the magnetoelastic energy density (Eq. ©). 
This is indicated by the relatively small quantitative dif- 
ferences in the calculated magnetoelastic coefficients for 
the various cation arrangements. On the other hand, 
the magnetoelastic coefficients of CFO show a stronger 
dependence on the specific cation arrangement and the 
orientation of the applied strain, so that the cubic ap- 
proximations is less justified in that case. In addition, 
the overall magnetoelastic response is much stronger in 
CFO than in NFO. 

Both of these observations can be understood from the 
d 7 electron configuration of the Co 2+ cation, which leads 
to stronger spin-orbit effects compared with the d 8 con- 
figuration of Ni 2+ . In the latter, the orbital magnetic 
moment is strongly quenched by the dominant octahe- 
dral component of the crystal-field, and the system is less 
sensitive to additional crystal field components of lower 



symmetry. In contrast, the orbital moment is not fully 
quenched by the octahedral crystal field for the d 7 con- 
figuration of Co 2+ , and additional splittings, which are 
created by the different arrangements of the surrounding 
B site cations, can have much stronger effects on the elec- 
tronic ground state within the partially filled minority- 
spin tig orbital manifold. 

Both sign and magnitude of the calculated magne- 
tostriction constants agree well with available experimen- 
tal data. Even for CFO, where the calculated magne- 
tostriction depends more strongly on the specific cation 
distribution than for NFO, the resulting uncertainty is 
within the spread of available experimental data. 

Further experimental data for single crystals is there- 
fore required for a more accurate comparison. We note 
that a number of obstacles can in principle affect an accu- 
rate comparison between theory and experiment. Apart 
from potential influences of varying sample stoichiome- 
try, degree of inversion, and measuring temperature, the 
preparation of an ideal demagnetized state with an es- 
sentially random orientation of magnetic domains is rel- 
atively hard to achieve. For example, a state with 50 % 
of domains oriented parallel and 50 % of domains ori- 
ented antiparallel with respect to a certain axis would 
have zero magnetization but the magnetostrictive strain 
would already be saturated along that direction. Further- 
more, for systems with very strong magnetic anisotropy, 
such as e.g. CFO, it can be very difficult to achieve full 
saturation along the hard direction^ Other sources of 
disagreement between theory and experiment could be 
due to the neglect of higher order terms in the energy 
expression ([2]), 2 or most likely due to deficiencies in the 
exchange correlation potential used in the DFT calcu- 
lations. However, based on the currently available ex- 
perimental data it can be concluded that the GGA+Z7 
method used in the present work is sufficiently accurate 
for further investigation on the effects of cation distribu- 
tion, degree of inversion, and stoichiometry on the mag- 
netostrictive properties of spinel ferrites. 

Our work thus provides a sound basis for future in- 
vestigations of magnetostriction and anisotropy in spinel 
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ferrites as well as for future first principles studies of mag- 
neto-electric coupling in artificial multifcrroic hctcrostruc- 
tures containing either CFO or NFO in combination with 
ferroelectric and/or piezoelectric materials. 
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